Constellation Plots: Neocortex, 2nd Trimester

All Samples, All cell types

Load data

---

• Neocortex_v3.RData Not necessary if just plotting on already calculated PCA, UMAP, clusters.

# Full dataset:
load("../data/cfc4b2_neocortex_v3.RData")

# 40k random subset (toy dataset)
# ncx.40k <- read_rds("constellation_plots/ncx_v3_40k.rds")

# 100k random subset of excitatory neuron lineage only.
ncx.exn.100k <- read_rds("../data/exn_lineage/toy/ncx.v3.exn.sub_100k.rds")

• Cluster / marker tables.

ncx.clusters <- read_delim("../tbls/83d19_Neocortex_allindividuals_combinedclusters_v1.txt",
                              "\t", escape_double = FALSE, trim_ws = TRUE)

ncx.markers <- read_delim("../tbls/8f0fc_Neocortex_subset1_clustermarkers_combo2.txt", 
                  "\t", escape_double = FALSE, trim_ws = TRUE)

---

scrattch.hicat steps to build constellation plots.

1. Define cells to build the plot from.

[All cells: 404.2K]

a. All clusters / cell types:

Keep only cells with combo2 $combined.cluster.2 annotation.

b. Excitatory neuron lineage only:

Keep only cells with combo2 $combined.cluster.2 annotation AND whose $combined.cluster.2 annotation belongs to excitatory neuronal lineage classes.

b.1. Cells in exn linage 100k cell toy object:

2. Build dataframes for constellation plots.

1. cl

2. rd.dat

3. cl.df

4. rd.cl.center Find cluster centers from UMAP coordinates

5. cl.center.df

Join cl.df and rd.cl.center into cl.center.df for input into get_KNN_graph.

6. knn.cl

Attaching package: ‘matrixStats’

The following object is masked from ‘package:dplyr’:

    count

3. Make constellation plot

---

2. Find DE genes between connected clusters.

---

Neuron_3 cells with nearest neighbors in Neuron_8 cell.cl.counts is a matrix of all cells and their counts of nearest neighbors in each cluster.

Find cells in cluster_a with bearest neighbors in cluster_b

nn.counts <- x[[cluster_b_col]]
length(x$cell.name)

[1] 10663

sum(nn.counts > 0)

[1] 6103

sum(nn.counts == 0)

[1] 4560

(median.nnCounts <- median(nn.counts[nn.counts > 0]))

[1] 5

Distribution of neighbor counts in cluster_b for cells in a given cluster_a. Use this to find the point at which cells will be split into comparison groups.

Compare cells above and below the median count of cluster_b neighbors.

cells <- list()
cells$above.median <- x %>% filter(get(cluster_b_col) > median) %>% pull(cell.name)
cells$below.median <- x %>% filter(get(cluster_b_col) < median)  %>% pull(cell.name)

s.obj <- SetIdent(s.obj, cells = cells$above.median, value = 'nn_ct_above_med') %>% 
          SetIdent(cells = cells$below.median, value = 'nn_ct_below_med')

table(s.obj@active.ident) %>% as.data.frame.matrix()
markers <- list()
markers$nn_above_med <- FindMarkers(s.obj, slot = "scale.data", 
                                  # cells.1 = cells$above.median, cells.2 = cells$below.median,
                                  ident.1 = "nn_above_med", ident.2 = "nn_below_med", 
                                  logfc.threshold = 0)

Calculate enrichment ratio, filter genes w. adj p-value < 0.05, sort table. Positive values indicate that the feature is more highly expressed in the first group.

markers.tmp <- markers$nn_above_med %>% rownames_to_column("gene") %>%
  mutate(enrich.ratio = pct.1 / pct.2,
         gene.score = avg_diff * enrich.ratio,
         across(.cols = where(is.numeric), .fns = round, digits = 5)
  ) %>%
  filter(p_val_adj <= 0.05) %>% 
  select(gene, pct.1, pct.2, enrich.ratio, avg_diff, gene.score, p_val_adj) %>%
  arrange(desc(gene.score))

write_tsv(markers.tmp, "../out/DEgenes_neuron3_vs_neuron8.tsv")

Make heatmap:

Extra functions